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Abstract 



' We describe a Fortran 90 program to calculate population ratios of atomic levels. 

The program solves the equations of statistical equilibrium considering all possible 
bound-bound processes: spontaneous, collisional or radiation induced (the later ei- 
^T) ' ther directly or by fluorescence). There is no limit on the number of levels or in 

. the number of processes that may be taken into account. The program may find 

a wide range of applicability in astronomical problems, such as interpreting fine- 
structure absorption lines or collisionally excited emission lines (such as coronal 
■ emission lines) in the spectra of several objects, and also in calculating the cooling 

, rates due to collisional excitation. 
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PROGRAM SUMMARY 

Title of program: PopRatio 
Catalogue identifier: 

Program obtainable from: |ht t p : / / www . iagusp . usp . br / ~ alexsilv /p opr at io 



Computer for which the program is designed and others on which it has been 
tested: 

Computer: Intel Pentium Pro 200 MHz PC 
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Installation: Institute Astronomico e Geofisico, Universidade de Sao Paulo. 

Operating system under which the program has been tested: DOS, Windows 
95. 

Programming language used: Fortran 90. 

Memory required to execute with typical data: 
2,055 kbytes 

No. of bits in a word: 32 

No. of processors used: One 

Has the code been vectorised or parallelized?: No 

No. of bytes in distributed program, including test data, etc.: 
981 058 

Distribution format: 
Keywords: 

Level populations, statistical equilibrium equations, fine structure lines, coUi- 
sionally excited lines, coronal lines, cooling rates. 

Nature of physical problem: 

The purpose of this program it to calculate population ratios of atomic levels. 
The excited level population ratios may be used to infer the physical conditions 
from observed fine-structure absorption lines or coUisionally excited emission 
lines (such as coronal emission lines). Another possible use is in calculating 
coohng rates due to coUisional excitation of low-lying levels. 

Method of Solution: 

The necessary atomic data is read from a separate input file. Next, the rates 
for all the bound-bound processes involved are evaluated and the system of 
equilibrium equations is solved. Several processes may be taken into account: 
spontaneous, collisions with an arbitrary number of particles, excitation or 
stimulated emission induced by radiation fields (either directly or by fluores- 
cence). Built in radiation fields provided are: a black body (such as the cosmic 
microwave background radiation), the UV radiation field of the Galaxy, the 
UV background of all QSOs and the hot halo model radiation field. Moreover, 
an arbitrary user-defined radiation field may also be included. 

Restrictions on the complexity of the problem: 

None. The code can handle an arbitrary number of processes and levels. The 
required memory is dynamically allocated in run time. 
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Typical running time: 

The following table gives the running times (in seconds) of the testcases pro- 
vided in section 4: 



System library Distributed libraries 



test case ^1 


0.06 


0.06 


test case #2 


0.06 


0.10 


testcase #3 


0.06 


0.06 


testcase #4 


0.05 


0.07 


testcase #5 


0.12 


0.14 



The listed values correspond to averages over 10 executions on a PC Pentium 
Pro 200 MHz. The column labelled " System library" gives the running times in 
case the routines from the MSIMSL library are used. The figures given in the 
next column correspond to the distributed version of PopRatio, that includes 
routines from the BLAS, LAPACK and PPPACK hbraries. The testcases were 
run using unoptimized BLAS routines. 

LONG WRITE-UP 



1 Introduction 

In most astrophysical environments, common atoms or ions are found pri- 
marily in their ground states, with only a negligible fraction in excited levels. 
However, if some excitation mechanism is present, then a small or even sig- 
nificant fraction of atoms/ions will be found in their excited states. If the 
excitation mechanism is not so intense to keep the atoms/ions in LTE condi- 
tions, then the amount of atoms/ions found in excited states may be used to 
infer the physical conditions in their environments. 

The fraction of excited atoms/ions may be infered from the observed spec- 
trum properties. For example, the ratios of excited level populations may be 
deduced from the column density ratios of fine structure absorption lines seen 
in absorption clouds towards QSOs or in diffuse clouds towards brigth stars in 
the Galaxy [1]. They may also be infered from intensity ratios of coUisionally 
excited emission lines (such as coronal emission lines). 

However, in order to trace the diagnosis curve giving the fraction of ex- 
cited atoms/ions as a function of the physical conditions in the medium, one 
must usually account for several excitation mechanisms [2]. Moreover, some 
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atoms/ions may have a complex enough electronic structure to require the 
inclusion of several levels in the calculation, such as iron. 

This paper describes a code to calculate the population ratios of excited lev- 
els of a given atom or ion accounting for an arbitrary number of levels and 
excitation mechanisms. 

Our code may find a wide range of applicability in astronomical problems, 
such as infering the physical conditions from fine-structure absorption line 
and coUisionally excited emission line diagnosis. Another possible use is in the 
calculation of cooling rates due to collisional excitation. We provide several 
examples in the testcases described bellow. 



2 Method of Solution 

2. 1 The equations of statistical equilibrium 

In order to calculate the level populations of a given atom or ion, we make 
two basic assumptions : 

(1) The rates of processes involving ionization stages other than the atom 
or ion being considered (such as direct photoionization or recombination, 
charge exchange reactions, collisional ionization, etc.) are slow compared 
to bound-bound rates. 

(2) All transitions considered are optically thin. 

To calculate the population of some level z, we must take into account all 
possible processes that will (de)populate it: 

(^fi ■ 

= ^ (processes that populate) — ^ (processes that depopulate) ,(1) 

where Ui is the volume density of atoms or ions in level i. 

Therefore, in steady state regime the sum over all processes that populate 
level i will be balanced by the sum over all processes that depopulate level i. 

Assuming that the two conditions listed above are met, this can be written 
(see, for instance, Rybicki and Lightman [3]): 

^ rij [Aji + BjiUji + ^ n'q^A = ^ [Aj + B^jUij + ^ n\1j J , (2) 

j \ k J j \ k / 
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where we have considered all possible bound-bound processes, i.e., sponta- 
neous, radiation-induced and collisionally-induced. The lefthand side of eq. 
(2) is the sum over all processes that populate level i from the other levels j, 
whereas the righthand side is the sum over all processes that depopulate level 
i to levels j. 

Aij is the transition probability of spontaneous decay from level i to level j. 
For i < j, Aij = 0. 

Bij are Einstein coefficients, related to the transition probabilities by: 



Bii ■■ 



Snh (Ei - EjY 



R — 

ji — ij ) 

9j 



(3) 



for i > j, and Ba = 0; h is Plank's constant, Ei is the energy of level i 
(expressed in cm"-*^) and Qi is the statistical weight of level i. 

Uij is the spectral energy density of the radiation field integrated along the 
line profile (pi, of the transition from level i to level j : 

Uij = Uji = / u^(f)^du ~ / {v - Vij) du = (uij) , (4) 



with Uii = 0; Uij is the frequency of the transition and we have assumed that 
the radiation field does not vary significantly along the line profile. 

In eq. (2) we have also considered the effect of collisions; is the volume 
density of the particle inducing the transition, the main collision partners 
usually being k — e~ , p'^ , , He^ , H2, ... , depending whether the medium is 
primarily ionized or neutral. 

q^j is the collision rate for the transition from level i to level j induced by the 
collision partner k . These coefficients are the cross-sections for the related 
process aij convolved with a Maxwellian distribution of velocities / (v) , mak- 
ing these quantities suitable for astrophysical applications (see, for instance, 
Osterbrock [4]): 



/I 00 
1 8kT f 
va {v) f (v) dv = j <^ii (e) ee'^de, 



(5) 



for the deexcitation rates {i > j); k is Boltzmann's constant, T is the kinetic 
temperature, fj, is the reduced mass of the system and e is the collision partner's 
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kinetic energy. 

Excitation and deexcitation rates are related by tiie principle of detailed bal- 
ance: 

Qij = — e qji, (6) 

9i 



with — 0. When the interaction is coulombian, as in collisions with elec- 
trons, it is convenient to express the cross-section in terms of the collision 
stregth fly, defined by: 



e - 7 

Snme gi 



where m is the electron's mass. Substituting this in eq. (5) yields: 



I 1 ^ , / e A 8.629 10-^ 7i,- 



with T expressed in K and 7.^ is defined by the integral in cq. (8) and is called 
Maxwellian- averaged collision stregth. Typically ^ij is a slowly varying function 
of T, of order unity. However, for neutral atoms it may vary for several orders 
of magnitude. 

These arc the basic parameters needed to solve eq. (2). If we consider our 
model ion to be composed of n levels, then we must solve a linear system of 
n — 1 equations in order to calculate the relative population ratios. 



2.2 Fluorescence 



Sometimes we have the situation in which the atom or ion is in one of its first 
m lower-lying levels and then it is photoexcited to some higher-lying level A 
(e.g., by an UV radiation field). Next it decays - either spontaneously or by 
stimulated emission - back to some different level among the first m levels. 
We call this process fluorescence. 

In this case it is easy to eliminate the level A from the linear system of equations 
(2), reducing its order by one. 

Initially, for notation purposes let us define: 

Kj^j = BijUij (9) 
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The statistical equilibrium equation (2) for some level I belonging to the first 
m lower-lying levels is: 

nx {Axi + Kxi) + ---^ niKix + • • • , (10) 



where we have written only the terms involving the level A. Whereas the 
equation for level A will be: 

m rn 

^ mKix = nxY. {Axi + Kxi) . (11) 
1=1 1=1 

Solving for in eq. (11), substituting in eq. (10) and introducing the indirect 
excitation rate from level i to level j through level A as: 

= K ^Aj + Kxj , , 



with r<^j = 0; we can now rewrite eq. (10) as, 



$:n,rj, + --- = n,^rf^. + ---. (13) 

3=1 J=l 



And we have eliminated the equation involving nx- Extending this reasoning 
to eliminate a whole set of higher-lying levels is straightforward. One simply 
replaces the indirect excitation rates in eq. (13) by the corresponding total 
indirect excitation rates: 

r^.^Er.^- (14) 

2.3 Solving the system of statistical equilibrium equations 

Let us now proceed to set the system of equations to calculate the relative 
population ratios. 

For notation purposes, let us define: 

Qij = Aij + Kij + + E ^""4 ' (15) 

k 

SO that eq. (2) reads (taking into account the possible effect of fluorescence): 

^rijQji^ Hi^Qij, (16) 



7 



or equivalently in matrix notation: 



Qij 

j 

Qn i^^j) (17) 

Hi 

0. 

After rewriting the system to have the solution expressed in terms of the 
populations relative to the ground level rii we are left with: 



Qw (18) 

Qj+i,i+i, ioriy^j; 
-EjQi+i,j: fori = j. 

Eq. (18) is the linear system of equations that is solved by PopRatio. 



3 The package PopRatio 

3. 1 The structure of PopRatio 

PopRatio is divided in four modules: PRECISION, PHYSMATH, FIELDIN- 
TENSITIES and POPRATIO. Next, we describe each one in greater detail. 

3.1.1 PRECISION 

This module defines single/double precision kind parameters and sets the ac- 
curacy to which the fioating point operations are to be done throughout the 
entire package. The user may change the precision of floating point operations 
by simply editing the WP parameter: 

integer, pcirameter :: WP = DP ! working precision definition 

The default value is set to correspond to the obsolescent Fortran 77's double 
precision data type. 



P-R^O 
P — 

P = 
Ri = 



M-X 



X; 



ni 

h — —Pi+1,1 
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3.1.2 PHYSMATH 



This module contains definitions of tt and some pliysical constants, as well as 
some mathematical routines. 

The mathematical routines are: 

• SWAP: subroutine that swaps any two integer numbers; 

• POL: function that evaluates a polynomial at a given point, 

• LINTERPOL: function that interpolates linearly a set of points. 

3.1.3 FIELDINTENSITIES 

This module contains routines that give the spectral energy densities (in 
erg/cm^/Hz) of some particular radiation fields taken from the hterature. The 
routines take as input the wavelength in A at which the field is to be evaluated 
(as well as some additional parameters required by each particular model) and 
interpolates linearly the table of points using routine LINTERPOL. 

The radiation fields currently available are: 

• FGAL: this function returns the radiation field of the Galaxy, as given by 
Gondhalekar et al. [5] ; 

• FUVB: this function returns the UV background radiation field of all QSOs, 
as given by Madau et al. [6]. It also takes as input the redshift at which 
the radiation field is to be evaluated Z and a flag MODEL. M0DEL=1 
corresponds to a cosmological model with = 0, M0DEL=2 to one with 
go = 0.5 and M0DEL=3 corresponds to go — 0-5 as revised by [6], 

• FHOTHALO: this function returns the radiation fleld of the hot halo model 
proposed by Viegas and Friaga [7] to explain the source of ionization of 
Lyman Limit QSO absorbtion hue systems. It takes as input two flags: R 
and T. R is the distance from the center of the galaxy at which the radiation 
flcld is to be evaluated, R=l corresponds to 10 kpc, R=2 to 30 kpc and R=3 
to 100 kpc. T is the age of the galaxy, T=l corresponds to 0.206 Gyr and 
T=2 to 0.3644 Gyr. 

3.1.4 POPRATIO 

This is the main module that actually computes the rate coefficients for the 
processes described in section 2 (after reading the related atomic data from 
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Table 1 

Correspondence between global variables in module POPRATIO and mathematical 
quantities. 



Variable Mathematical quantity 



E(i) 








A(i i) 




Bfi 






Uij 


CD(k)%q(i,j) 




gam(i,j) 


r 


matrix(i,j) 


Mij 


indep(i) 


h 


X(i) 





a separate input file) and then solves the system of equations of statistical 
equilibrium. 

It begins defining a few useful constants to the evaluation of the rate coeffi- 
cients and three new types: 

• TABLE: this type is a table of points, Y vs. X; 

• TCOLLISIONALDATA: this type encloses all necessary information per- 
taining a certain coUisional process, such as the collision partner name, the 
collision rates, etc., 

• THLEVEL: this type encloses all necessary information pertaining some 
of the ji higher- lying levels described in section 2.2, such as the energy, 
statistical weight, transition probabilities, radiative rates, etc. 

This module also declares several global variables that are visible to the main 
program. Table 1 shows the correspondence between most important variables 
and the mathematical quantities given in section 2. 

We also have the following routines: 

• SKIP_COMMENTS: subroutine used by INITPOPRATIO to skip the com- 
ments from the atomic data input file (see section 3.2 bellow); 

• INITPOPRATIO: subroutine that reads in the atomic data values stored in 
the input file. It takes as input a unit number for the input file (defined via 
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an open statement in the main program), the number of levels NBLEVELS 
to be taken into account in the calculation, and wavelength bounds to the 
radiation field intensities, LAMBDAI and LAMBDAF. The last two argu- 
ments are optional parameters, and if supplied by the user INITPOPRATIO 
will not load in any upper level connected to the ground term by multiplets 
whose wavelengths lie outside of this range (the default values are 0.0_WP 
and huge(1.0_WP)). INITPOPRATIO does not load into RAM memory all 
the data stored in the input file, only the necessary space to take into ac- 
count the desired number of levels is allocated in run time. (See also section 
3.2 below); 

• FINISHPOPRATIO: subroutine that deallocates the dynamical variables 
allocated by INITPOPRATIO; 

• POPULATIONRATIO: this is the key subroutine that sets the rates for 
all the processes and calculates the population ratios. It takes as input the 
kinetic temperature of the gas T (in K), the redshift Z, a vector N containing 
the volume densities (in cm~^) of the collision partners entered in the atomic 
data input file and two optional parameters: BETA and F. 

The redshift is used to set the temperature of the Cosmic Microwave 
Background Radiation (CMBR), as given by the following relation predicted 
by the standard Big Bang cosmology (see, for instance, Kolb and Turner 
[8]): 

T^Toil + z), (19) 

where Tq = 2.728±0.002 K is the current value of the CMBR temperature 
[9,10]. The effect of the CMBR blackbody radiation field is then included 
in the energy densities: 

Uij = p{i^ij,T) = — — ^ , (20) 

where c is the speed of light. Alternatively, there are models that gener- 
alize the standard temperature law of the CMBR (19) by introducing a free 
parameter /? [11]: 

r = To(l + z)'-^ (21) 

where /3 is within the range < < 1. If this parameter is not entered via 
the optional parameter BETA, it is assumed zero. POPULATIONRATIO 

calls function TCMBR to evaluate the temperature of the CMBR, and next 
function PLANCK to evaluate its contribution to the energy densities. If 
the user does not wish to take the CMBR into account, the redshift Z should 
be entered with the value -1.0_WP. 
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POPULATIONRATIO also adds to the energy densities the contribution 
from a radiation field defined by the user in function URAD. The optional 
parameter F multiplies the radiation field defined in URAD by a constant 
factor /. If F is not entered, it is assumed 1. 

LTEDEV: this function returns the module of the largest relative deviation 
of the population of the lowest NL levels from the population expected from 
LTE conditions characterized by a temperature T_LTE (in K); 

URAD: this function returns the spectral energy density of a user-defined 
radiation field. The user may create his/her own subroutine and then call 
it in the body of this function. 

PLANCK: this function returns the spectral energy density of a blackbody 
radiation field (20). It first checks whether the argument of the exponential 
does not exceed the natural logarithm of the largest number in the mashine 
representation in order to avoid overflow (the spectral energy density is set 
to zero in case this happens). 

TCMBR: function that returns the temperature of the CMBR (21); 

RATE_EXTRAPOL: this function is called by POPULATIONRATIO when- 
ever the kinetic temperature T lies outside the range allowed by some table 
of points of coUisional rates entered in the atomic data input file (see section 
3.2 below). The user must supply a subroutine to extrapolate the table of 
points and modify this function to call it, otherwise a warning message will 
be issued telling which particular table of points needs to be extrapolated. 
An analytical fitting to some transition may be implemented by entering an 
unphysical temperature range in the input file, so that whatever value of T 
is entered, POPULATIONRATIO wiU always call this function that may be 
used to call another function containing the fitting (such as CI_GAM and 
OLPROTON, see example in section 3.2 below); 

CI_GAM: function that returns the MaxweUian averaged collision strengths 
for transitions - induced by collisions with electrons - involving the lowest 
five levels of C° (based upon analytical fittings taken from the literature 
[12,13]), 
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• OLPROTON: function that returns the deexcitation rates by coUisions with 
protons for transitions involving the lowest three levels of 0° (based upon 
analytical fittings taken from the literature [14]). 



3.1.5 Other packages used by PopRatio 

PopRatio uses several routines from the BLAS [15], LAPACK [16] and PP- 
PACK [17] libraries. The routines needed are distributed together with PopRa- 
tio in case the user does not have them installed in his/her system. 

We use the Fortran 90 interface to the LAPACK routines described by Don- 
garra et al. [18]. To call the PPPACK routines we have written a new module - 
SPLINE - containing Fortran 90 interfaces to the routines and a new function 
that calls them: SPLINE3. The PPPACK routines were shghtly modificated 
in order to support both real and double precision arguments. 



3.2 The atomic data input file 

The user must supply all necessary atomic data used by PopRatio in a sepa- 
rate input file, that will be read by INITPOPRATIO. We provide input files 
containing atomic data pertaining five atoms/ions of interest in fine structure 
absorption line diagnosis of QSO absorbers: C°, C"*", Si"*", 0° and Fe"*". Although 
in these files we list all the references in the literature where the atomic data 
were taken, we do not discuss the particular models adopted here, since this 
will be the subject of a forthcoming paper (Silva and Viegas, in preparation). 

Next we describe the general syntax of the input file. As an example we take 
the file OLdat, that contains the atomic data of neutral oxygen. 

The atomic data information is distributed along several blocks, each one 
starting with the "7^" character. The user is free to insert any comments 
before this symbol. The first block: 

# 
01 

is the name of the atom/ion being considered. This string is stored in the 
global variable SPECIES_NAME. It may be up to 10 characters long. The 
next blocks of data are: 

# 
3 
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this is the number of levels that belong to the ground term. In atomic oxy- 
gen these levels are ^P2,i.o- PopRatio will take into account fluorescence for 
transitions among these levels. This value is stored in the global variable 
NGROUND. 



# 
5 

this is the maximum number of levels that may be taken into account. The 

user must make sure to enter enough data to allow PopRatio to consider this 
number of levels. If INITPOPRATIO is called with a value for NBLEVELS 
higher than this, a warning message will be issued. 



# 



0.0 


5 


1 


2s2 2p4 3P2 


158.265 


3 


2 


2s2 2p4 3P1 


226.977 


1 


3 


2s2 2p4 3P0 


15867.862 


5 


4 


2s2 2p4 1D2 


33792.583 


1 


5 


2s2 2p4 ISO 



these are the values for the energies and statistical weights for the levels. The 
first column is the energy relative to the ground level (in cm~^), and the 
second is the statistical weight for the level. The last two columns numbering 

and naming the levels are not read by INITPOPRATIO, they are just entered 
for user reference. The energies must be sorted in ascending order. The energies 
and statistical weights are stored in global variables E(i) and G(i), respectively. 



# 
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1 


2 


8.865E-5 


1 


3 


1.275E-10 


2 


3 


1.772E-5 


1 


4 


6.535E-3 


2 


4 


2.111E-3 


3 


4 


6.388E-7 


1 


5 


2.945E-4 


2 


5 


7.909E-2 


4 


5 


1.124E0 



This block contains the transition probabilities. The number at the beginning 
is the number of transitions listed. The first column is the lower level index j for 
the transition, next are the higher level index i and the transition probability 
Aij (in s~^). The transition probabilities are stored in the global variable A(i,j). 
As the transition probabilities are read, INITPOPRATIO automatically sets 
the Einstein coefficients by eq. (3) and stores them in the global variable 
B(i,j). There is no need to enter null transition probabilites, since transition 
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probabilities not entered are assumed to be zero. 



135 








76794.798 


3 


1 


3.30E+08 


76794.798 


3 


2 


1.97E+08 


76794.798 


3 


3 


6.54E+07 


96225.049 


3 


1 


9.31E+07 


96225.049 


3 


2 


5.56E+07 


96225.049 


3 


3 


1.85E+07 


97488.378 


3 


1 


1.99E+06 


97488.378 


3 


2 


2.97E+07 


97488.378 


3 


3 


3.95E+07 


97488.448 


5 


1 


1.79E+07 


97488.448 


5 


2 


5.35E+07 



These are the transitions involving the /i higher- lying levels of section 2.2 and 
the NGROUND ground term levels. The number at the beginning is the num- 
ber of transitions listed. The first column is the energy relative to the ground 
level (in cm'^) of the higher-lying level, the second is its statistical weight, 
the third is the lower-lying level index (<NGROUND) for the transition and 
the fourth is the transition probability (in s~^). The energy levels must be 
entered in ascending order. If some higher-lying level was already loaded in by 
INITPOPRATIO in the fourth block (and therefore will be explicitly included 
in the system of statistical equilibrium equations), then it will not be taken 
into account when POPULATIONRATIO calculates the indirect excitation 
rates, GAM(i,j). However its contribution will be added in TOTGAM(i,j), in 
case the user wishes to assess it. 

# 
4 

This is the number of collision partners. This value is stored in the global 
variable NPARTNERS and is the dimension of the coUisional data type vector 
CD(i). 

Next we have one block for each one of the coUisional partners: 

# 

electron 

5 

10 

3 7 

50. 100. 200. 500. 1000. 2000. 3000. 

2 1 8.34E-04 1.26E-03 1.79E-03 2.58E-03 3.35E-03 4.61E-03 5.92E-03 



15 



3 1 3.23E-04 5.00E-04 7.34E-04 1.12E-03 1.49E-03 2.08E-03 2.64E-03 
3 2 2.74E-07 8.92E-07 2.78E-06 1.05E-05 2.38E-05 4.62&05 6.98E-05 



The first line is the coUisional partner name, with a maximum of 10 characters 
long. This string is stored in the global variable CD(i)%PARTNER_NAME. 
The second line is the index of the highest-lying level that may be pop- 
ulated by collisions with this partner; it is stored in CD(i)%MAXLEVEL. 
The third line is the total number of transitions hsted below; it is stored in 
CD(i)%NINTERPOL. For each transition we have a table of points contain- 
ing the collision rates (or Maxwcllian-averagcd collision strcgths, if the partner 
name is 'ELECTRON') as a function of the kinetic temperature. PopRatio in- 
terpolates this table of points by a cubic-spline in log-log space (linear-log, 
in case the partner name is 'ELECTRON'). PopRatio also uses eq. (6) to 
calculate the inverse process rate. 

In the next line we have the partial number of transitions that follows and the 
number of data points entered in the tables. The next line lists the temperature 
values (in K) for the transitions that follows. Then we have the collision rates 
(in cm^s~^, or Maxwellian-averaged coUision strengths), the first two values 
entered giving the i and j indexes for the transitions {i — > j). If the partner 
name is electron then the order of these indexes are irrelevant, otherwise i must 
come first. This is repeated until the sum of the partial number of transitions 
equals the total number of transitions listed. This is to allow greater flexibility, 
since rates for the various transitions might be taken from different sources 
in the literature and therefore have been calculated in different temperature 



values. 






# 






proton 






3 






3 






3 2 








1.E30 


1.E30 


2 1 


l.E-30 


l.E-30 


3 1 


l.E-30 


l.E-30 


3 2 


l.E-30 


l.E-30 



This is an example of how the user can implement an analytical fitting. The 
transitions induced by collisions with protons are entered, but with an un- 
physical temperature range. Since the temperature passed to POPULATION- 
RATIO will always lie outside this range, it will call RATE_EXTRAPOL that 
is then used to call OLPROTON, containing the analytical fittings. 

The last two blocks in Ol.dat give the collision rates for transitions induced 
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by collisions with neutral hydrogen and neutral helium. 



3.3 Installing and running PopRatio 



3.3.1 Setting the accuracy of floating point operations 

The user must set the accuracy of the floating point calculations by editing 
the WP kind parameter in module PRECISION (see section 3.1.1). 

There are only two options available, corresponding to the single/double pre- 
cision Fortran 77's data types. This restriction applies because PopRatio uses 
F77 routines from the BIAS, LAPACK and PPPACK libraries. The Fortran 
90 interfaces to these routines will automatically call the correct routines suit- 
able to the selected precision. 

The default precision is set to correspond to the double precision Fortran 
77 data type. 



3.3.2 Setting the radiation field 

The user must supply the radiation field he/she is interested in by editing the 
body of function URAD in module POPRATIO. 

For example, in order to include the FGAL and FUVB radiation fields de- 
fined in module FIELDINTENSITIES the user should change the body of 
this function to: 

urad = Fgal(lambda) 
if (z/=-1.0_WP) then 

urad = urad -|- Fuvb(lambda,z,3) 
end if 

If the user does not intend to use any of the radiation fields defined in module 
FIELDINTENSITIES then it does not need to be included in the project, and 
the corresponding use Fieldlntensities statement may be deleted from module 
POPRATIO. 

The default is a null radiation field: 
urad = 0.0_WP ; 
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3.3.3 Writing the main program 

A minimal main program should include an open statement specifying the 
atomic data input file, a call InitPopRatio statement loading in a model for the 
atom/ion being considered and a call FinishPopRatio statement to deallocate 
the dynamical variables, freeing up RAM memory. 

Typically between the calls to INITPOPRATIO and FINISHPOPRATIO we 
will have several calls to POPULATIONRATIO under different physical con- 
ditions. After each call to POPULATIONRATIO the rates for all the processes 
involved and the population ratios will be acessible via the global variables 
given in table 1. 

The calculation will proceed faster in case the physical conditions are not very 
different from the physical conditions in last call to POPULATIONRATIO. 
If the same value for the kinetic temperature T is passed in two successive 
calls to POPULATIONRATIO, then the collision rates CD(i)%q(i,j) will not 
need to be re-computed. Similarly, if we have the same values for Z, BETA 
and F parameters, the radiation field energy densities U(i,j) and the indirect 
excitation rates GAM(i,j) will remain unchanged. 

Finally, a note of caution: if the user is interested in the collision rate Qij by a 
given coUisional partner, then he/she should not call POPULATIONRATIO 
with its density set to zero, since the collision rates will not be calculated in 
this case. 

In the next section wc give some examples of main programs to tackle common 
applications. 



4 Applications and examples of using PopRatio 

4.-1 Testcase #1 

In this example we calculate the population ratios of the ^Pj ground levels 
of atomic carbon as a function of neutral hydrogen density. These population 
ratios are a useful density estimator [19] and are needed to interpret ratios 
of UV C I fine structure absorption lines observed in damped Lyman-a QSO 
absorbers [20-22] , for example. 

We take into account collisions with several collisional partners - electron, 
proton, neutral hydrogen, para/ortho molecular hydrogen, neutral helium -, 
tying their densities to the neutral hydrogen density n^o according to: 
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rie. 



Up — O 



no-H2 = '^''^HO 



(22) 



We also take into account the effect of CMBR at a redshift z = 1 and the 
radiation field of the Galaxy. 

In the high density hmit, function LTEDEV is called to show that the popu- 
lation ratios differ from LTE conditions by less than 4% . 

To run this testcase the user should include the radiation field of the Galaxy 
in the body of function URAD in module POPRATIO: 

urad = Fgal(lambda) 



4.2 Testcase #2 



Now we want to determine the neutral hydrogen density from the population 

ratio (deduced from the observed column density ratios of UV lines in the 
spectrum) assuming the same physical conditions as in the previous example. 

This example exploits the fact that the diagnosis curve traced in the first 
example is a monotonically increasing function of density. 

To run this testcase the user should include the radiation field of the Galaxy 
in the body of function URAD in module POPRATIO: 

urad — Fgal(lambda) 



4.3 Testcase #3 



This example reads in the neutral hydrogen density determined in the previous 
example and, assuming the same physical conditions as above, determines 
the rates for all the processes involved in order to access what are the most 
important excitation mechanisms. 

The results show that the most important excitation mechanism are collisions 
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by neutral hydrogen, with the CMBR making a significant contribution for 
the ^Pi level excitation. 

To run this testcase the user should include the radiation field of the Galaxy 
in the body of function URAD in module POPRATIO: 

urad = Fgal(lambda) 
4-4 Testcase #4 

We now illustrate how to use PopRatio to calculate coUisionally excited emis- 
sion line intensity ratios. 

As an example we take the UV intercombination multiplet 2s2p^ ^P — > 2s^2p ^P° 
of C II at 2325 A. Intensities ratios of lines belonging to this multiplet may 
serve as a useful indicator of electronic densities in the range 10^ < Ue < 
10^^ cm~^ [23,24]. These lines have been observed in the spectra of a variety 
of astronomical objects: planetary nebulae [25], giant stars [26,27], symbiotic 
stars [28] and in the solar chromosphere and transition region [29] . 

The emissivity of a line is given by [4] : 

e (Xij) = Hi Aij hi^ij . (23) 

We calculate the following emissivity ratios as a function of electronic density: 




We take into account collisions by electrons and protons and fluorescence 
induced by a black body radiation field of temperature = 4000K attenuated 
by a geometric dilution factor / = 0.5 (which might be representative of a 
stellar chromosphere). 

In the high density limit function LTEDEV is called to show that the level 
populations differ from LTE conditions by less than 0.4%. 
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To run this tcstcase the user must include the black body radiation field in 
the body of function URAD in module POPRATIO: 

urad = Planck(4.0E3_WP,c/lambda*1.0E8_WP) 



4-5 Testcase #5 



In this example we illustrate how PopRatio may be used to calculate the 
coohng rate due to coUisional excitation of low-lying levels of a given atom or 
ion. 

As an example we take the ion Fe"*", because iron is an astrophysically abundant 
element. The ion Fe"*" may be the dominating ionization stage in low ionization 
regions, and may be an important gas coolant [30]. 

The cooling rate is given by [4] : 

^■max 

L{T) = J2J2^^ A, hly^J . (25) 

i=2 j<i 



Rewriting this in terms of the total Fe"*" density npe+ = J2i and the popula- 
tion ratios Xj = 

L(T) -1 

= {1 + Y.X,) Yin ^-1 (^^ - , (26) 

'^Fe+ i i=2 j<i 



with the energies Ei expressed in cm ^. The Fe+ density will depend on its 
fractional abundance and on the iron elemental abundance: 

riFe+ = riH (27) 

riFe riH 



Here we calculate the cooling function, defined as the right hand side of eq. 
(26). 

Since Fe"*" has a comphcated electronic structure, several levels must be taken 
into account in the calculation. We employ a 16-level model ion, allowing us 
to calculate the coohng function for electronic densities as high as 10^ cm~^. 

To run this testcase the user does not need to modify function URAD, since 
no fiuorescence transitions are loaded in. 
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5 Concluding remcirks 



We have described a program to calculate atomic level population ratios. The 
code may find a wide range of applicability in astronomical problems, as we 
have illustrated in the testcases. 
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